## setup

# clean environment
rm(list=ls())

# load packages
library(ri2)

# load data
dat <- read.csv2("data_experiment_outcomes_csv.csv",
                encoding = 'UTF-8')

# outcome variable to numeric
dat$response <- as.numeric(as.character(dat$response))

# filter out untreated
dat <- dat[!is.na(dat$response),]


## randomization inference

# ri: immigrant background treatment effect
effect.perm<-NULL
set.seed(1)
for(i.perm in 1:50000){
  # Permute the treatment status
  treat.perm<-sample(dat$immback, size = 1522, replace = F)
  # Calculate the treatment effect
  effect.perm[i.perm]<-mean(dat$response[treat.perm==1])-mean(dat$response[treat.perm==0])
}
sum(abs(effect.perm) >= abs(mean(dat$response[dat$immback==1]) - mean(dat$response[dat$immback==0]))) / 50000

# ri: immigrant background, rightist parties
effect.perm<-NULL
set.seed(1)
for(i.perm in 1:50000){
  # Permute the treatment status
  treat.perm<-sample(dat$immback, size = 1522, replace = F)
  # Calculate the treatment effect
  effect.perm[i.perm]<-mean(dat$response[treat.perm==1 & dat$left==0])-mean(dat$response[treat.perm==0 & dat$left==0])
}
sum(abs(effect.perm) >= abs(mean(dat$response[dat$immback==1 & dat$left==0]) - mean(dat$response[dat$immback==0 & dat$left==0]))) / 50000

# ri: immigrant background, leftist parties
effect.perm<-NULL
set.seed(1)
for(i.perm in 1:50000){
  # Permute the treatment status
  treat.perm<-sample(dat$immback, size = 1522, replace = F)
  # Calculate the treatment effect
  effect.perm[i.perm]<-mean(dat$response[treat.perm==1 & dat$left==1])-mean(dat$response[treat.perm==0 & dat$left==1])
}
sum(abs(effect.perm) >= abs(mean(dat$response[dat$immback==1 & dat$left==1]) - mean(dat$response[dat$immback==0 & dat$left==1]))) / 50000

# ri: immigrant background, conditional on sender being a partisan
effect.perm<-NULL
set.seed(1)
for(i.perm in 1:50000){
  # Permute the treatment status
  treat.perm<-sample(dat$immback, size = 1522, replace = F)
  # Calculate the treatment effect
  effect.perm[i.perm]<-mean(dat$response[treat.perm==1 & dat$partisan==1])-mean(dat$response[treat.perm==0 & dat$partisan==1])
}
sum(abs(effect.perm) >= abs(mean(dat$response[dat$immback==1 & dat$partisan==1]) - mean(dat$response[dat$immback==0 & dat$partisan==1]))) / 50000

## ri: partisan
effect.perm<-NULL
set.seed(1)
for(i.perm in 1:50000){
  # Permute the treatment status
  treat.perm<-sample(dat$partisan, size = 1522, replace = F)
  # Calculate the treatment effect
  effect.perm[i.perm]<-mean(dat$response[treat.perm==1])-mean(dat$response[treat.perm==0])
}
sum(abs(effect.perm) >= abs(mean(dat$response[dat$partisan==1]) - mean(dat$response[dat$partisan==0]))) / 50000

## ri: interaction
effect.perm<-NULL
set.seed(1)
for(i.perm in 1:50000){
  # Permute the treatment status
  treat.perm<-sample(dat$immback, size = 1522, replace = F)
  treat.perm2<-sample(dat$partisan, size = 1522, replace = F)
  # Calculate the treatment effect
  effect.perm[i.perm]<-(mean(dat$response[treat.perm==1 & treat.perm2==0])-mean(dat$response[treat.perm==0 & treat.perm2==0])) - (mean(dat$response[treat.perm==1 & treat.perm2==1])-mean(dat$response[treat.perm==0 & treat.perm2==1]))
}
sum(abs(effect.perm) >= abs((mean(dat$response[dat$immback==1 & dat$partisan==0]) - mean(dat$response[dat$immback==0 & dat$partisan==0]))-
                              (mean(dat$response[dat$immback==1 & dat$partisan==1]) - mean(dat$response[dat$immback==0 & dat$partisan==1])))) / 50000

